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Abstract. We numerically calculate the evolution of second order cosmological perturbations 
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1 Introduction 



The arrival of more and better data in the last couple of years has fundamentally changed 
theoretical cosmology. Until only a few years ago it was sufficient to use linear order the- 
ory to analyse the data and calculate first order observables, such as the power spectrum; 
higher order theory wasn't required. The new data, in particular the Cosmic Microwave 
Background (CMB) anisotropy maps as provided by the Wmap and the Planck satellites 
[1, 2], make it possible to go beyond linear order and derive higher order observables such 
as the bispectrum from the data. Considerable effort has therefore been spent on extracting 
observable signatures beyond the power spectrum and to extend cosmological perturbation 
theory beyond linear order. 

So far calculations of higher order observables such as the bispectrum, or its popular 
parametrisation /nl,j have relied heavily on approximations. This meant usually assuming 
large scales, that is large compared to the "horizon", and imposing slow-roll [3]. However, 
these restrictions not only constrain the validity of the results, but also limit the number of 
models that can be studied, excluding many interesting cases. For example, even in linear 
theory it is already well known that interesting and observable effects occur when slow-roll 
is violated at the end of inflation [4] . 

Second order perturbations play a crucial role in our quest to understand the non-linear 
physics of the early universe. Previous works by us and collaborators [5, 6] have derived the 
equations of motion for second order perturbations in the long wavelength approximation, 
the slow-roll approximation and also the full non-slow- roll case for all scales [7]. 

Previously in Ref. [8] we numerically solved the evolution of second order perturbations 
under the assumption that the source term, which consists of convolutions of first order 
perturbations and their derivatives, is calculated in the slow-roll approximation. We showed 
that the source term closely follows the form of the first order power spectrum. In this paper 
we go beyond the slow roll approximation to calculate the full source term equations for 
a single scalar field. This full treatment is needed in cases where slow-roll is broken. The 
results of the updated calculation are consistent with those of our previous work for the ^m 2 ip 2 
potential for which the slow-roll approximation is an exceptionally good one. To go beyond 
slow-roll in this paper, we study "step" and "bump" potentials, for which slow-roll is broken 
as 1 77v | > 1 around the feature (rjy denoting one of the slow-roll parameters). Step potentials 
have been widely used to create features in the power spectrum of inflationary perturbations 
which might more accurately match the observed power spectrum from WMAP [9-14]. In 
this paper we will follow the potential forms described by Chen et al. [15, 16]. In addition 
to being able to calculate the second order field perturbation in this case, the result of the 
full equations seems to be smoother at very early times when the perturbations are highly 
oscillatory (due to a larger damping term, which is lost in slow-roll). 

The final goal of this continuing work is a numerical calculation of the curvature per- 
turbation at second order for all length scales. This will probe effects both inside and outside 
the horizon in a way that is not possible using other methods, for example the 5N formalism 
[17-21], which a priori uses large scales only. The applications of such a calculation are 
many and varied and range from an investigation of the non-Gaussian nature of inflationary 
perturbations to the exploration of higher order effects such as the sourcing of gravitational 
waves at second order [22] and vorticity generation [23]. 

In the current paper we have only considered the effect of a single scalar field. One 
goal of future work is to extend the numerical code to deal with more than one field. In 
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a multifield system the curvature perturbation is not the only important observable with 
isocurvature components also playing a significant role. Calculating the field evolution is the 
only way to gain access to these isocurvature results, in a way that would not be possible 
if we had concentrated purely on the curvature perturbations. With this goal in mind we 
have fashioned the numerical calculation in terms of a scalar field instead of the curvature 
perturbation which might otherwise have been considered a more observationally relevant 
quantity in the single field case. We have developed the numerical code considerably from 
that used in Ref. [8] . The change from the slow-roll to full equations adds significant complex- 
ity to the calculation but this has been countered by optimising the time sensitive functions 
and implementing more of them in a compiled language. The numerical code has been re- 
leased under an open source license (Modified BSD license) and is available for download [24] . 

In Section 2 we present the Klein-Gordon equation for second order cosmological per- 
turbations. This is the central part of the code described in Section 3. The results we have 
obtained are in Section 4 and discussion of these and future goals is contained in Section 5. 
Throughout the text we use 8irG and label conformal time derivatives with a prime. We 
assume a flat Friedmann-Robertson- Walker background throughout this paper. 



2 Second Order Perturbations 
2.1 Klein Gordon Equations 

Cosmological perturbations beyond linear order have been studied in depth in recent times 
and a brief summary can be found in Appendix A or at length in Refs. [25, 26] (see Ref. [27] 
for an example of earlier work in this area). In particular the Klein-Gordon equations for 
the background, first order and second order perturbations are given in Eqs. (A.15)-(A.17). 
These equations are given in real space, however, and in general the dynamics of the scalar 
field become clearer in Fourier space. Following Refs. [3] and [7] we will write 5(f(k l ) for the 
Fourier component of 5<p(x l ) such that 



6<p{v, x l ) = J & s k5^{k l ) exp(ikix l ) , (2.1) 

where k l is the comoving wavenumber. In Fourier space, the closed form of the first order 
Klein-Gordon equation (A. 16) then transforms into 



tipple)" + 2n5ip l {k i y + k^s^k* 



+ a 2 



5<pi(k i ) = 0, (2.2) 



where ipo is the background homogeneous scalar field, 6<pi is the first order scalar field 
perturbation, a is the FRW scale factor, T~L = a' /a is the conformal Hubble parameter and 
V is the potential of the scalar field with derivatives w.r.t ip denoted by etc. The second 
order equation (A. 17) requires a careful consideration of terms that are quadratic in the first 
order perturbation. Terms at second order of the form ((5c/?i(x 1 )) 2 require the use of the 
convolution theorem (see for example Ref. [28]). We use convolutions of the form 

/(^(a-i) JL J dW^'-y-o'l/fp'lfltf). (2-3) 
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The full closed form, second order Klein-Gordon equation in Fourier space is then given by 
[7] 

<y 2 '(^) + 2H5<f' 2 {k i ) + k 2 5i P2 (k i ) 

.2 
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(2.4) 



The source term S contains all the products of 6<pi in real space which require convolution 
integrals. Terms which contain gradients of 8(pi include additional factors of k and q. The 
form of S is given by [7] 
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where we have defined the parameter Q = a 2 (8irGVoip' /'H + V (p ) for convenience. Calculating 
Eq. (2.5) is the main task of the numerical simulation described in Section 3. 



2.2 Equations for Numerical Calculation 

In the previous section the governing equations of the second order system were given in 
terms of conformal time. A more appropriate time variable for the numerical simulation is 
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the number of e-foldings, Af = log(a/a; n it). This is measured from ai n ; t , the value of a at the 
beginning of the simulation. We take the value of the scale factor today, do, as ao = 1 and 
calculate a at the end of inflation by connecting it with ao (see discussion in Ref. [29] or for 
example Eq. (3.19) in Ref. [3]). We also assume that reheating was instantaneous at the end 
of inflation. We will use a dagger (') to denote differentiation with respect to Af. Derivatives 
with respect to conformal time, r\, and coordinate time, t, are then given by 
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(2.6) 
(2.7) 



dt dr] dAf 

respectively, where H = dlna/dt and % = aH. The background and first order equations, 
written in terms of the new time variable Af, are 



0. 



(2.« 



and 
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6< Pl (k i ) = 0. (2.9) 



The three-dimensional convolution integral J dq 3 can be rewritten in spherical coordinates 
q,8,u; where q = \q l \. There is no u dependence in the source term integral but \k l — q l \ 
and the factors of h % and q % are dependent on 9 and this dependence is made explicit below. 
Eqs. (2.4) and (2.5) must be written in terms of Af, with the 9 dependent terms grouped 
together, in order to set up the numerical system completely at second order. In the slow-roll 
case there were only four different 9 dependent terms, here labelled A-T> following Ref. [8]: 
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(2.10) 

The non-slow-roll source term in Eq. (2.13) that we are now considering requires the use of 
three further 9 integrals in addition to those in Eq. (2.10), which are 
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It is worth noting that the term sm s (6)/\k l — q l \ 2 tends to zero in the limit where k = q and 
# — t- 0. The second order Klein- Gordon equation in e- folding time is 
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where the full source equation is given by 
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For comparison, the slow-roll expression for the source term is [8] 
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The increased length of calculation in going from the slow-roll source term in Eq. (2.14) to the 
full version in Eq. (2.13) can be clearly seen. The numerical complexity is not significantly 
greater, however, once the three terms in Eq. (2.11) have been calculated. In the next section 
the implementation of the numerical scheme to calculate the source term in Eq. (2.13) and 
solve the second order Klein-Gordon equation is discussed. 



3 Code Implementation 

We described in Ref. [8] a numerical system which calculated the slow-roll source term as 
given in Eq. (2.14). We have improved the system and the full source term calculation has 
been implemented. This code has now been released under an open source license [24]. The 
basic structure of the system is still the same as outlined in Ref. [8] and follows the form 
laid down by Ref. [30-32]. Firstly the background fields are evolved to pinpoint the end of 
inflation when en = —H^/H = 1. The first order Klein-Gordon equation (2.9) is solved 
using a fourth order Runge-Kutta scheme with the initial conditions specified by the Bunch- 
Davies vacuum Eq. (3.1). 1 Following this first order stage the source term is calculated at all 
necessary timesteps. 2 In the last stage of the numerical calculation, the source term result is 
used to evolve the second order perturbation modes using Eq. (2.12). 

The major advance reported in this paper is the use of the full (non-slow-roll) source 
equation in the third step above. With this advance, seven different terms need to be calcu- 
lated at each time step. These vary in k, q and 9 and integrals over q and 9 are performed, 
again at each time step, for each k value. In order to perform these integrals we implement 
cutoffs at large and small values of k. These cutoffs and the integration calculations in general 
are described in detail in Refs. [8, 29] and we do not discuss them further here. The seven 
different terms are functional forms of 9 given in Eqs. (2.10) and (2.11). As the first order 
results are tabulated at specific k values, it is necessary to interpolate the results for k val- 
ues between these points. The increased complexity of the source equation clearly increases 
the number of calculations at each time step and hence the overall execution time. The 
original slow-roll calculation was numerically intensive so the effect of almost doubling the 

1 The parameters used in the potential have previously been fixed by fitting the first order power spectrum 
to the best fit WMAP value of Viz — 2.45 x 10 -9 , five e-foldings after horizon crossing. 

2 This stage is naively parallelisable as the calculation at each timestep is independent. To shorten execution 
time the source term can be calculated for selected k modes only instead of the full range. Note that the first 
order calculation still needs to be run with a large range of k modes in order that the convolution integral can 
be performed. 
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number of calculations could have been dramatic. In fact, optimisation of the code and the 
translation of some parts into compiled modules has made the increase in complexity man- 
ageable. By selecting only particular modes of interest, for example the WMAP pivot scale 
fewMAP = 0.002MPC- 1 = 5.25 x l(r 60 Mp L , the source term calculation can be significantly 
shortened. 3 

In order to solve the ODEs given in Eqs. (2.8), (2.9) and (2.12) initial values of 
ipo, Scpi, 5if2 and their first derivatives must be specified. For the background field initial val- 
ues of (/?o and ip^ are chosen so that the inflationary period lasts approximately 60 e-foldings 
after the scales of interest in the CMB exit the horizon. In the case of the quadratic potential 
this requires choosing super-Planckian initial field values, for example </?o,init = 18Mpl- 

The first order perturbation initial conditions are chosen by assuming that sufficiently 
far inside the horizon the perturbation modes are in the Bunch-Davies vacuum state. To 
implement this early time condition in the numerical system we follow Salopek et al. [30] 
who initialise the first order modes at a time when k/aH for the mode equals some arbitrary 
factor. In keeping with Ref. [30] we choose this factor to be 50. The first order initial values 
are then calculated as 



where r] is the conformal time again. 

The situation of the second order initial conditions is different. At the initial times 
when the Bunch-Davies conditions are suitable the perturbations are expected to be highly 
Gaussian. In this paper we are interested in the production of second order effects by the 
evolution of the first order modes and we make no assumptions about the existence of second 
order perturbations before the simulation begins. Therefore, we set the initial condition for 
each second order perturbation mode to be 5<p2 = 0, and 8<p\ = at the time when the 
corresponding first order perturbation is initialised. 

A numerical solution for the second order perturbation equation will contain a homo- 
geneous solution and a particular solution. As stated above we have chosen the initial values 
for the second order field to be zero. On their own these initial conditions do not remove 
this homogeneous solution from the result for 5if2 in general. In order to do this, and keep 
only the particular solution to the equation, it is necessary to ensure that the homogeneous 
solution is the trivial (0, 0) solution throughout the evolution. 

In order to only report the particular solution of the second order differential equation 
Eq. (2.12) we have added a ramping term to the source term S. This ramp interpolates 
between and 1 for about an e-folding of time around the time of the initialisation of the 
modes. By starting the solution of Sf2, <5v4 a ^ 0, and setting the source term to at this 
time through the use of the ramp, the solution for S(p2 will consist only of the inhomogeneous 
part. The homogeneous solution of Eq. (2.12) i.e., the solution without the source term 
present, is the same form as the solution for 5ip± as can been seen by comparing Eqs. (2.12) 
and (2.9). Figure 1 shows the effect that incorporating the ramp has on the source term at 
early times for the fewMAP mode. The value of 5^2 (&wmap) is initialised to be equal to zero 

3 The run time of the source term calculation for a single k mode is approximately two hours on a relatively 
modern CPU. Using a cluster of computing nodes shortens the naively parallelisable source term calculation 
to a further degree. 
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(a) The ramp function used to remove the homoge-(b) The addition of the ramp term to S changes its 
neous solution for 8<p2, here shown around the initial- value at the initialisation time from the original value 
isation time for the scale fcwMAP- (blue solid line) to the ramped value (green dashed 

line), here shown for scale fcwMAP. 

Figure 1: The ramping function used with the source term results. 



at 64.3 e-foldings before the end of inflation, when the ramp value is still zero. The ramp 
is added to the source term for each mode at the respective initialisation time and all the 
results below were generated using a ramped source term. 

4 Results 

4.1 Comparison with Slow-Roll Results 

The full second order code has been run with a number of different potentials. Firstly in 
order to check the consistency of the full equation, the standard slow-roll quadratic potential 
was used. Figure 2a shows that the results of the full system are very similar to the slow-roll 
results, as expected for this potential. The additional terms in the source equation Eq. (2.13) 
subdue some of the oscillatory noise evident in the slow-roll solution at early times when the 
mode is inside the horizon. 

The second order result is similar in both cases, however there is an appreciable increase 
in the amplitude of the second order modes when the full equations are used. This is likely to 
be a result of the reduced oscillations mentioned above. The second order values are plotted 
in Figure 2b. 

The results for background and first order perturbations are robust under small devia- 
tions in the initial conditions or the mass parameter in the potential. However, the source 
term and second order results are not so robust, with any slight variation in the initial value 
of ipo or the mass parameter translating to sizable changes in the magnitude of the results. 
For example, a small change in the initial field value of Ay?o,init = O.OIMpl or in the mass of 
the inflaton of Am = 1 x 1CP u Mpl leads to large differences in the source term and second 
order results and Figure 3 shows the effects of these changes. 

4.2 Step and Bump Potentials 

Beyond the standard quadratic model, a more interesting potential to consider is one with a 
feature at a particular scalar field value [9]. Following Chen et al. [15, 16] we have used both 
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(a) The source term results at all time steps (b) The second order results around the feature 



Figure 2: A comparison of the results from the full equations (red solid line) and the slow-roll 
source term (green dashed line) for scale &wmap- 




Figure 3: The source term and second order results for models with slightly perturbed mass 
(green dashed line) and cp^ (blue dotted line) compared to the standard quadratic model 
(red solid line) for scale &wmap- 



a step and a bump potential. The step potential is a modified \rr?ip 2 potential of the form 



2 2 

-m (p 



1 + c I tanh 



d 



(4.1) 



and the bump potential is given by 
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where the parameters c, d, ip s and ip^ control the height, width and central point of the feature 
respectively. The step potential used here has an extra (—1) term compared to the one used in 



Ref. [15]. This extra term ensures that V s — > hm 2 (p 2 at early times instead of beginning with 
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(a) The potential V(<p) for the step, bump and stan-(b) The slow-roll parameter r\y for the step, bump 
dard quadratic potentials around the feature. and quadratic potentials. 

Figure 4: The potential and r\y slow-roll parameter for the step (red solid line), bump (green 
dashed line) and quadratic (blue dotted line) potentials. 
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0.0018 


0.0005 
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0.022 


0.01 




14.84 M PL 


14.84 M PL 



Table 1: The values of the parameters in the step and bump potentials. 



a greater amplitude than the standard quadratic potential and only equalling the quadratic 
value exactly at the feature point. The evolution of the Hubble parameter H is different 
for both the bump and step potentials compared to the quadratic potential. This changes 
the value of a at the end of inflation and therefore the value a; n it used at the start of the 
run. In order to compare models with similar initial conditions the following results are for 
runs where ai n it has been made equal in each case. In physical terms this amounts to small 
redefinitions of the value of a today away from unity, or the consideration of slightly different 
physical scales today. Due to the strong dependence on initial conditions as shown above, 
comparison of numerical results is better facilitated by this fixing of ai n ;t than comparing 
models with different values of ai n it . 

Figure 4a shows the step and bump potentials at the relevant <p values. The features are 
quite shallow for this choice of parameters but can be tuned to be stronger. The values used 
in the code are given in Table 1 . With these values the slow- roll approximation is temporarily 
violated as \rj\ gets large around the feature as shown in Figure 4b. We have run the step 
potential model with a "full" step, corresponding to c = 0.0018, a "half" step with c = 0.0009 
and finally with c set equal to zero. As a "sanity check" the c = run gives back the same 
results as the standard quadratic potential with no feature. The bump potential model has 
also been run with a "full" bump for which c = 0.0005, a "half" bump with c = 0.00025 and 
a "zero" bump. 

The first order spectrum for the three cases for the step potential is plotted in Figure 5 
and the cases for the bump potential in Figure 6. At first order the effect of the feature in both 
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(a) The first order results at all time steps 



(b) The first order results around the feature 



Figure 5: The first order results for the full (red solid line), half (green dashed line) and 
zero (blue dotted line) step potentials for scale fcwMAP- 
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(b) The first order results around the feature 



Figure 6: The first order results for the full (red solid line), half (green dashed line) and 
zero (blue dotted line) bump potentials for scale &wmap- 



potentials is localised around a particular ip value and therefore J\f value. The half step and 
bump deviations are smaller than the full ones, however the ratio of the amplitude changes 
are not symmetric around the central point of the feature. In the source term calculation 
the presence of a feature makes a great difference in the result. Figures 7 and 8 compare the 
results for the step and bump potentials again with a full, half and zero feature. For the step 
potential the magnitude of the source term deviates from the standard quadratic result quite 
far in advance of the feature. Interestingly the change in magnitude around the feature is 
almost equal in the full and half cases and is certainly not proportional to the parameter c in 
the way the first order results are. The magnitude of the source term is reduced compared to 
the quadratic result and this reduction continues beyond the feature. This reduction occurs 
before the feature when the step potential should be well approximated by the quadratic one. 
The higher order derivatives will be different however and in particular V^ w for the step 
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(a) The source term results at all time steps (b) The source term results around the feature 

Figure 7: The source term results for the full (red solid line), half (green dashed line) and 
zero (blue dotted line) step potentials for scale /cwmap- 
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where cp s is the value of the field at the feature. Unlike the third derivative of the quadratic 
potential this is not zero everywhere. The effect on the source term can be understood by 
examining Eq. (2.13) and observing that the first term is proportional to V^ w . In addition, 
the second term in the equation is proportional to U w and will not be constant as in the 
quadratic case. In contrast, the results for the bump potential are not affected beyond the 
small region around the feature, as shown in Figure 8. Again the change in magnitude does 
not seem to be proportional to the parameter c. Before and after the feature the result is 
indistinguishable from the quadratic case even though higher order derivatives of the potential 
are different to the values for the quadratic potential. 

At second order the differences in the models make themselves even more apparent. 
Figure 9 shows the second order solution for fcwMAP again for the full, half and zero step 
potentials. The amplitude of the step has a marked effect on the amplitude of the second 
order modes. The differences in the source term before the feature are carried over to the 
second order result. The magnitude of the second order result is much lower for the full step 
potential than the the quadratic result. When the amplitude of the step is halved the change 
in the magnitude is also reduced. The difference between the full and quadratic results is at 
least an order of magnitude and the detailed cause of this difference is a subject for further 
investigation. The results for the bump potential are shown in Figure 10. Here the effect of 
the feature is only evident around the bump, again carrying over the result from the source 
term values. 

4.3 Sub- and Super-Horizon Features 

For different k modes the feature occurs when the mode is either inside or outside the horizon. 
When the mode has already crossed the horizon the result is as in Figure 5. Modes which 
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(a) The source term results at all time steps (b) The source term results around the feature 

Figure 8: The source term results for the full (red solid line), half (green dashed line) and 
zero (blue dotted line) bump potentials for scale /cwmap- 
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Figure 9: The second order results for the full (red solid line), half (green dashed line) and 
zero (blue dotted line) step potentials for scale /cwmap- 



are inside the horizon when they encounter the feature in the potential are not affected in 
the same way. Of course the location of the bump can also be changed by varying the ip^ 
parameter in the potential. In Figure 11 the first order power spectrum and source term 
results are plotted for a potential where the bump feature is located inside the horizon 
and (/9b = 15.5Mpl and compared with the standard quadratic potential and the normal 
bump potential where </?b = 14.8Mpl. When the bump is sub-horizon (the red solid line in 
Figure 11) the first order results are slightly perturbed from the standard result, reaching a 
slightly altered magnitude after horizon crossing. This does not happen for the super-horizon 
bump (green dashed line) which asymptotes back to the quadratic result after the feature. 

The source term results also differ depending on the location of the bump. When the 
bump is sub-horizon the change in the magnitude of the source term is much more suppressed 
compared to when the bump is outside the horizon. In the second case the oscillations of 
nearly all the modes have already been dampened by the time the feature is encountered. 
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(a) The second order results at all time steps (b) The second order results around the feature 

Figure 10: The second order results for the full (red solid line), half (green dashed line) and 
zero (blue dotted line) bump potentials for scale &wmap- 



The convolution integral over the modes is then much more affected by the change in the first 
order perturbations around the feature. In contrast when the bump is sub-horizon, at least 
for the &WMAP scale, most of the other scales considered in the integral are still oscillating 
and the net effect is a small change in the magnitude of the source term. At second order 
the results are similar to first order. The sub- horizon bump slightly changes the magnitude 
away from the quadratic result, in this case reducing it. This change is kept beyond the hori- 
zon in contrast to the super-horizon case where the result asymptotes smoothly back to the 
quadratic one beyond the feature. Figure lid shows the ratio of the sub- and super- horizon 
bump second order results to the standard quadratic result. Following the sub- horizon bump 
the difference between the results is of the order of 2% whereas for the super-horizon bump 
the results are indistinguishable from the quadratic results after the feature region. 

In this section we have outlined the main results of our new numerical calculation. We 
have shown that the source term results using the full non-slow-roll equation (2.13) are more 
damped than the corresponding slow-roll results. We have demonstrated the new code by 
using feature potentials with a step and a bump added to the standard quadratic potential. 
Depending on the position of the feature and the form of the addition to the potential the 
effects of the feature can be seen in the source term and second order results beyond the 
feature itself. Whether the feature is sub- or super-horizon for a particular mode also affects 
the subsequent evolution. 



5 Discussion and Conclusion 



We describe in this paper the numerical solution of the full Klein-Gordon equation for a single 
scalar field at second order in cosmological perturbation theory. We use gauge-invariant 
variables in the flat gauge without imposing the slow-roll approximation or using the large 
scale limit. This is an extension of previous work, that relied on the slow-roll approximation 
to calculate the source term [8]. We have made this extended code publicly available (it can 
be downloaded from the website [24]). To validate and test the code, we studied several single 
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(c) The second order results around the feature, (d) The ratio of the second order results for the sub- 
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Figure 11: The first order and source term results for the bump potential when the bump 
is inside the horizon when ip^ = 15.5Mpl (red solid line) or outside the horizon when ip^ = 
14.8AfpL (green dashed line). The standard quadratic potential results are also shown as the 
blue dotted line and all the results are for the scale fcwMAP- 



field potentials, including step and bump extensions of the simple chaotic inflation potential, 
which violate one of the slow-roll conditions. 

We have shown that feature potentials affect the source and second order results and at 
least for the step potential these effects are apparent throughout the evolution of the modes 
and not just in a narrow region around the feature. This is consistent with the higher order 
derivatives of the potentials being non-zero, or at least having different values when computed 
without a feature. Is this a problem for feature potentials? At one level the feature potentials 
we have used are theoretical toy models which are not motivated by high energy physics. At 
first order in perturbation theory the tanh and cosh additions to the quadratic potential 
work to localise the deviation to a specific region. This is clearly not the case at second 
order for the step potential. To replicate this behaviour at second order, a restricted version 
of the step or bump potentials could be used. However, one main benefit of using tanh or 
cosh is that they smoothly asymptote to the original potential away from the feature. Any 
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attempted patching together of featureless potential and feature in a specific region would 
have to be carefully constructed to avoid discontinuities in the derivatives of the potential. 
In contrast the results for the bump potential are only affected in a narrow region around 
the feature. It is interesting to note, that a step-potential of the form (4.1) has recently been 
used to study the generation of magnetic fields in the early universe [33]. 

Somewhat surprising was the effect the small sub-horizon bump had on the evolution of 
the second order field fluctuation 5 if 2 on large scales, see Figure 11. Sub- horizon effects like 
this can only be studied by codes like the one described here, but would not be calculable by 
super- horizon codes, such as the recently proposed one in Refs. [34, 35]. However, it is not 
clear yet whether this change in bipi would translate into observational consequences, such 
as a change in /nl- We are cautiously optimistic that it might have an observable effect, but 
postpone a detailed calculation to future work. In this future work we plan to study more 
complicated and also more interesting potentials than in Section 4, and also hope to extend 
the code to allow the evolution of more than one field. In addition, to our knowledge there 
is currently no expression that exists in closed form for £2 or 7^2 which is valid on all scales 
without imposing slow-roll. Constructing such an expression should be possible and would 
allow for direct calculation of the evolution of £2 and we look forward to making progress 
on this issue. A recent work discussing the numerical calculation of /nl used lattice field 
theory simulations [36] . Whether in future this method will prove more efficient than the one 
discussed here will depend on detailed comparative analyses which are beyond the scope of 
the present article. Other recent work by Takamizu et al. [37] matched a perturbative solution 
for the curvature of a single field model inside the horizon to a non-linear solution outside the 
horizon which uses a gradient expansion. The authors claim this goes beyond the limitations 
of the 6N formalism and can also deal with temporary violations of slow-roll. It would be 
interesting in future work to compare the results of our approach with those of Takamizu 
et al. especially in the region where the matching of the solutions takes place. Gong, Noh 
and Hwang have also looked at higher order perturbations [38]. Their work concentrates on 
evolving the convolution kernels of the curvature perturbation but only for large scales with 
slow roll parameters equated to zero. The "Generalized Slow Roll" approximation can also 
be used to investigate models which break slow roll transiently around a feature [10, 39, 40]. 
Adshead et al. [40] have demonstrated calculations of the bispectrum of a single field model 
and observe good agreement with analytic results. However, this approach is limited to 
superhorizon scales only and can only be applied when there is a single degree of freedom. 
Another natural application and extension of our code is to apply it to other problems that 
have evolution equations of a similar form. As pointed out in the introduction, two immediate 
applications might be the numerical study of the generation of gravitational waves at second 
order [22] (see Ref. [41] for earlier numerical work) and the generation of vorticity [23]. In 
both cases the source term of the second order quantities is given by a convolution integral 
quadratic in first order quantities. 

Finally, the numerical calculation of the convolution integral requires cutoffs to be im- 
plemented at both large and small scales. We have explained the cutoffs we use explicitly 
in Refs. [8, 29] and refer the interested reader to the discussions contained therein. It will 
be interesting to implement in future work different cut-off schemes, as outlined for example 
in Refs. [42, 43]. Our current choice of cut-offs is pragmatic, however there may be physical 
considerations for adopting a different cut-off scheme, which would impact on the higher 
order perturbations and therefore affect the non-linear observables. 
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A Appendix 



The methods adopted for the study of first order perturbations can be extended at second 
order to find gauge invariant quantities. Recall that scalar quantities such as the inflaton 
field, cp, can be split into an homogeneous background, cpo, and inhomogeneous perturbations. 
Up to second order (p becomes 



ip(rj, x^) = <f (v) + fyiiv, x% ) + ^^(f?, x l 



(A.l) 



The metric tensor g^ u must also be perturbed up to second order. Here we consider 
only the scalar metric perturbations [26] : 
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90i 



-a 2 (1 + 20i + <fo) , 
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B x + -B 2 
1 2 



gij = a 2 [(1 - 2V>i - ip 2 ) Sij + 2Ex,ij + E 2tij ] 



(A.2) 



where 8%j is the flat background metric, <pi and (f> 2 are the lapse functions at first and second 
order, ipi and ip 2 are the curvature perturbations, and B\, B 2 , E\ and E 2 are the scalar 
perturbations describing the shear. As well as the first order transformation vector, there is 
a second order transformation vector and they are both given by 



(A.3) 



where the spatial vector part of the transformation has been ignored. 

The transformation of a second order scalar quantity (such as S(p 2 ) is given by [26, 44]: 

5<fi2 = Sf2 + ^0"2 + «1 (^0 a l + { P'o a 'l + 2 <Vl) + ( 2 <Vl + f'o a l) ti Pi, 1 > ( A - 4 ) 

where a tilde (~) denotes a transformed quantity. The metric curvature perturbation trans- 
formation at first order is straightforward, ipi = t/j\ — Tiai, but at second order it becomes 
more complicated [25, 26]: 

i^ + i 

where Xij is given by 
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and Bu and Cuj are denned as 

B\i = Bi i , Cuj = — V>i<% + Ei^j . (A. 7) 

Working in the uniform curvature gauge, where spatial 3-hypersurfaces are flat, implies that 

4>i = ^ 2 = E\ = E 2 = . (A.8) 

These relations can be used at first and then at second order to define gauge invariant 
variables. The first order transformation variables in the flat gauge satisfy a\ = tpi/7i and 
fix = —E\. At second order, for the transformation of scalar quantities, as in Eq. (A. 4), we 
require only This is found by using Eq. (A. 5) to have the form 



02 



^2 1 

n m 



■i.i 



(A. 



where the first order gauge variables have been substituted into Xij. 

The Sasaki-Mukhanov variable, i.e., the field perturbation on uniform curvature hyper- 
surfaces [45, 46], is given at first order by 



S<Pi + f^Vi 



(A.10) 



At second order the Sasaki-Mukhanov variable becomes more complicated [5, 44]: 
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(A.ll) 



where £ifl a t = —E\^. From now on we will drop the tildes and only refer to variables calculated 
in the flat gauge. The potential of the scalar field can also be separated into homogeneous 
and perturbed sectors: 



V(<p) 
5Vi 
6V 2 



V + 8Vx + -5V 2 , 
V w 5(fl + V<p5ip2 ■ 



(A.12) 

(A.13) 
(A.14) 



The Klein-Gordon equations are found by requiring the perturbed energy-momentum tensor 



Tuu to obey the energy conservation equation V nT^ v 



0. 



background field, ipo, the Klein-Gordon equation is 

+ 2H(f/ + a 2 V 
The first order equation is 

Sip" + ZHScp'x + 2a 2 V^4>\ - V 2 5ipi - ip' V 2 Bi - <ff <t>i + a 2 V )Lpip 5<pi 
and the second order one is given by 



(see for example Ref. [5]). For the 

(A.15) 



0, 



(A.16) 



6(p'i + 2H5(p 2 - V 2 5ip 2 + aV w ^2 + «V w (^i) 2 + 2a 2 V, v fa - <p' (V 2 £ 2 + <f>' 2 ) 



+ Vo^Lfc^if + 2 { 2 Wo + « 2 Kv) BijkB* + 40i (a 2 V >ipcp S<Pi ~ V 2 <Vi) 
= 0, 



(A.17) 
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where all the variables are now in the flat gauge. 

In order to write the Klein-Gordon equations in closed form, the Einstein field equations 
are also required at first and second order. These are not reproduced here, but are presented 
for example in Section II B of Ref . [7] . 
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